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We analyze the spectrum of harmonic flow, Vnipr) for n = 0-5, in event-by-event hydrodynamic 
simulations of Pb-t-Pb collisions at the CERN Large Hadron Collider (-^Sjvw = 2.76 TeV) with 
principal component analysis (PC A). The PC A procedure finds two dominant contributions to the 
two-particle correlation function. The leading component is identified with the event plane Vnipr), 
while the subleading component is responsible for factorization breaking in hydrodynamics. For vo, 

Ui, and V 3 the subleading flow is a response to the radial excitation of the corresponding eccentricity. 

By contrast, for V 2 the subleading flow in peripheral collisions is dominated by the nonlinear mixing 
between the leading elliptic flow and radial flow fluctuations. In the V 2 case, the sub-sub-leading 
mode more closely reflects the response to the radial excitation of £ 2 . A consequence of this picture 
is that the elliptic flow fluctuations and factorization breaking change rapidly with centrality, and 
in central collisions (where the leading V 2 is small and nonlinear effects can be neglected) the 
subsub-leading mode becomes important. Radial flow fluctuations and nonlinear mixing also play a 
significant role in the factorization breaking of U4 and 115. We construct good geometric predictors 
for the orientation and magnitudes of the leading and subleading flows based on a linear response to 
the geometry, and a quadratic mixing between the leading principal components. Finally, we suggest 
a set of measurements involving three point correlations which can experimentally corroborate the 
nonlinear mixing of radial and elliptic flow and its important contribution to factorization breaking 
as a function of centrality. 


I. INTRODUCTION 


tionally parametrized by r 2 {pTi,PT 2 ) [H], 


Two-particle correlation measurements are of 
paramount importance in studying ultrarelativistic 
heavy ion collisions, and provide an extraordinarily 
stringent test for theoretical models. Indeed, the 
measured two-particle correlations exhibit elliptic, 
triangular, and higher harmonics flows, which can be 
used to constrain the transport properties of the quark 
gluon plasma (QGP) [1, 2]. The remarkable precision 
of the experimental data as a function of transverse 
momentum and pseudorapidity has led to new analyses 
of factorization breaking, nonlinear mixing, event shape 
selection, and forward-backward fluctuations [3-8]. 
In this paper we analyze the detailed structure of 
two-particle transverse momentum correlations by using 
event-by-event (boost-invariant) hydrodynamics and 
principal component analysis (PCA) [9, 10]. Specifically, 
we decompose the event-by-event harmonic flow Vnipr) 
into principal components and investigate the physical 
origin of each of these fluctuations. This paper extends 
our previous analysis [10] for triangular flow at the LHC 
(Pb-I-Pb at ^Sjvjv = 2.76 TeV) to the other harmonics, 
n = 0-5. In particular, we demonstrate the importance 
of radial flow fluctuations for subleading flows of higher 
harmonics. 

Taking the second harmonic for definiteness, the two- 
particle correlation matrix of momentum dependent el¬ 
liptic flows, C 2 {pti,PT 2 ) = {V 2 {pti)V 2 iPT 2 )) Is tradi¬ 
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If there is only one source of elliptic flow in the event 
[for example if in each event V 2 (pt) = fiPT)s 2 with 
£2 a complex eccentricity and fipr) a fixed real func¬ 
tion of pr] then the correlation matrix of elliptic flows 
C 2 iPTi,PT 2 ) factorizes into a product of functions, and 
the £2 parameter is unity. However, if there are multi¬ 
ple independent sources of elliptic flow in the event, then 
the correlation matrix does not factorize, and the £2 pa¬ 
rameter is less than unity [11]. The r 2 parameter has 
been extensively studied both experimentally [3, 12-14] 
and theoretically [10, 11, 15, 16]. In particular, in our 
prior work on triangular flow we showed that factoriza¬ 
tion breaking in event-by-event hydrodynamics arises be¬ 
cause the simulated triangular flow is predominantly the 
result of two statistically uncorrelated contributions—the 
linear response to £3 [17] and the linear response to the 
first radial excitation of £3 [10]. The goal of the cur¬ 
rent paper is to extend this understanding of factoriza¬ 
tion breaking to the other harmonics. This extension was 
surprisingly subtle due to the quadratic mixing between 
the leading and subleading harmonic flows. 

Experimentally, it is observed that factorization break¬ 
ing is largest for elliptic flow in central collisions (see in 
particular Fig. 28 of Ref. [13] and Fig. 1 of Ref. [3]). 
Indeed, the r 2 parameter decreases rather dramatically 
from mid-central to central collisions. This indicates that 
the relative importance of the various initial state fluctu¬ 
ations which drive elliptic flow are changing rapidly as a 
function of centrality. The current manuscript explains 
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the rapid centrality dependence of factorization breaking 
in V 2 as an interplay between the linear response to the 
fluctuating elliptic geometry, and the nonlinear mixing of 
the radial flow and average elliptic flow. This quadratic 
mixing is similar to the mixing between and V 2 ,V 3 
[18-21], and this picture can be confirmed experimentally 
by measuring specific three point correlations analogous 
to the three plane correlations measured in the ^ 5 ,^ 2 , ^3 
case [4, 5]. 

To understand the linear and nonlinear contributions 
quantitatively, we will break up the fluctuations in hy¬ 
drodynamics into their principal components, and ana¬ 
lyze the linear and nonlinear contributions of each prin¬ 
cipal component to the simulated harmonic spectrum. 
The sample of events and most of the PCA methods are 
the same as in our previous paper [10], and therefore 
in Sects. IIA and IIB we only briefly review the anal¬ 
ysis definitions, and the key features of simulations. In 
Sect. lie we discuss the strategy for constructing the 
best linear predictor for leading and subleading flows. 

The second part of our paper. Sect. Ill, contains in¬ 
dividual discussions for each harmonic flow. First, we 
discuss radial flow fluctuations in Sect. Ill A and then 
demonstrate their importance in generating subleading 
elliptic flow in Sect. IIIB 1. In Sect. IIIC we briefly de¬ 
scribe our PCA results for direct and triangular flows. 
Finally, in Sect. HID we discuss the quadrangular and 
pentagonal flows and how the nonlinear mixing of lower 
order principal components adds to these flows. We put 
forward some experimental observables in the discussion 
in Sect. IV. A catalog of figures showing the main results 
of PCA for each harmonic is given in the Appendix. 


II. PCA OF TWO-PARTICLE CORRELATIONS 
IN EVENT-BY-EVENT HYDRODYNAMICS 

A. Principal components 

PCA is a statistical technique for extracting the dom¬ 
inant components in fluctuating data. In the context 
of flow in heavy ion collisions it was first introduced in 
Ref. [9], and then applied to the analysis of triangular 
flow in our previous paper [10]. Here we review the es¬ 
sential definitions. 

The event-by-event single particle distribution is cus¬ 
tomarily expanded in a Fourier series 

7 /\r 

-^ = Vo{pT) + Y.^n{PT)e-^^‘^+R.c., ( 2 ) 

^ n—1 

where dp = dy dpT dp denotes the phase space, p is the 
azimuthal angle of the distribution, and H.c. denotes 
Hermitian conjugate. Vnipr) is a complex Fourier co¬ 
efficient recording the magnitude and orientation of the 
nth harmonic flow, without the typical normalization by 
multiplicity. 

PCA is done by expanding the covariance matrix of 
two-particle correlations (which is real, symmetric, and 
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FIG. 1. Factorization ratio ri(pT\,PT 2 ) [Eq. (7)] for trian¬ 
gular flow and its approximations with principal components 
(PCs) in central collisions (0-5%). 


positive semi-definite) into real orthogonal eigenvectors 
V^^\pt), 

Cn{pTl,PT2) = ) (iC.p,, - (R:.p,J )) 

= (3) 


where is the square root of the 

eigenvalue times a normalized eigenvector. The eigen¬ 
value records the squared variance of a given fluctuation. 

The principal components Vn^\pT), V^\pt), ■ ■ ■ of a 
given event ensemble can be used as optimal basis for 
event-by-event expansion of harmonic flow 

Rn(pT) - {Vn{pT)) = + ■■■■ 

( 4 ) 

The complex coefficients are the projections of har¬ 
monic flow onto principal component basis and record 
the orientation and event-by-event amplitude of their re¬ 
spective flows. Principal components are mutually un¬ 
correlated 

(5) 


Typically the eigenvalues of Cn{pTi,PT 2 ) are strongly 
ordered and only the first few terms in the expansion are 
significant. Often the large components have a definite 
physical interpretation. We define the scaled magnitude 
of the flow vector Vn‘^\pT) as 




/ {yn°'\pT)) dpT 

f {dN/dpr)'^ dpT ' 


( 6 ) 


which is a measure of the size of the fluctuation without 
trivial dependencies on the mean multiplicity in a given 
centrality class. 

The leading flow vector Vn^\pT) corresponds to fluc¬ 
tuations with the largest root-mean-square amplitude, 
while subsequent components maximize the variance 
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in the remaining orthogonal directions. This yields a 
very efficient description of the full covariance matrix 
CniPTi,PT 2 ) and factorization ratio 


rn{PTl,PT2) 


Cn(,PTl,PT2) 


■\/Cn (PTI , PTI ) Cn {PT2 , PT2 ) 


< 1 - ( 7 ) 


rniPTi,PT 2 ) is bounded by unity within hydrodynam¬ 
ics [11]. By truncating series expansion of the covariance 
matrix [Eq. (3)] at two or three principal components we 
can approximate Cn(pTiiPT 2 ) and rniPTi,PT 2 )- Trun¬ 
cating at the leading term would constitute complete flow 
factorization, i.e., r„ = 1. The factorization matrix for 
triangular flow is shown in Fig. 1. We see that at low mo¬ 
mentum pt < ‘2. GeV just two principal components are 
sufficient to describe momentum dependence of factoriza¬ 
tion ratio r^. Analogous decompositions of two-particle 
correlations into principal components exist for all har¬ 
monics and all centralities. Interpreting these large flow 
components physically is the goal of this paper. 


B. Simulations 


We used boost-invariant event-by-event viscous hydro¬ 
dynamics to simulate 5000 Pb-Pb collisions at the CERN 
Large Hadron Collider (LHC) = 2.76 TeV) in 

fourteen 5% centrality classes selected by impact param¬ 
eter. The initial conditions are based on the Phobos 
Glauber Monte Carlo [17] with a two-component model 
for the entropy distribution in the transverse plane. We 
used a lattice equation of state [22] and a “direct” pion 
freeze-out at = 140 MeV. The results presented here 
were simulated using a shear viscosity to entropy ratio 
oi r]/s = 0.08. Qualitatively similar results are obtained 
with rj/s = 0.16, with the most important differences 
discussed in Sect. HIB2. The same ensemble of events 
was used in Ref. [10], which provides further simulation 
details. 


C. Geometrical predictors 

We will construct several geometric predictors for the 
leading and subleading flows following strategy outlined 
in Ref. [20] . Keeping the discussion general, let Ci^pj-ed 
a geometric quantity which predicts the event-by-event 
amplitude and phase of the corresponding flow . For 
example, for the leading n = 3 component the triangular¬ 
ity £ 3,3 (defined below) is an excellent choice for ^ 3 pj.ed' 
The geometric predictors are designed to maximize the 
correlation between a particular flow signal and the ge¬ 
ometry. Specihcally, the predictors maximize the Pear¬ 
son correlation coefficient between the event-by-event 
magnitude and orientation of ath principal component. 


and the geometrical predictor Ci'^pred 


max = 


lJa)p*(a) \ 
Snpred/ 


//c(“)t*Ca)\/c(“) ^ 

Y XS.Ti Sn /\‘5*n pred‘5>n pred/ 


( 8 ) 


We constructed several predictors for the flow by 
assuming a linear relation between the flow and the geom¬ 
etry. The simplest predictor consists of linear combina¬ 
tions of the Hrst two eccentricities of the initial geometry. 
These are defined as 


£n,n 


£n,n+2 




i?” 


R 




(9a) 

(9b) 


where the square brackets [] denote an integral over the 
initial entropy density for a specific event, normalized by 
the average total entropy 5tot- Rrms = \J ([?'^]) is the 
event averaged root-mean-square radius. Note that our 
definitions of £n,n and £n,n +2 are chosen to make the 
event-by-event quantities £n,n and £n,n +2 linear in the 
fluctuations. In this notation, the geometric predictor 
based on these eccentricities is 


?i°pred — ^n,n + £n,n+2, (10) 

where ci is adjusted to maximize the correlation coeffi¬ 
cient in Eq. (8), and the overall normalization is irrele¬ 
vant. While the first two eccentricities provide an excel¬ 
lent predictor for the leading flow, they do not predict 
the subleading flow very well. This is in part because the 
radial weight is too strong at large r. 

More generally, one can dehne eccentricity as a func¬ 
tional of radial weight function p{r): 


n 

rms 


( 11 ) 


It is the goal of this paper to hnd the optimal radial 
weight function p(r) for predicting both leading and sub¬ 
leading flows. For the subleading modes p{r) will have 
a node, and thus e„{p(r)} will measure the magnitude 
and orientation of the first radial excitation of the geom¬ 
etry [10]. 

To find the optimal radial weight we expand p{r) in 
radial Fourier modes 


9"r7l 

P{r) = ^ Wb-^Ju{hr), ( 12 ) 

where Jn{x) is a Bessel function of order n, Wb are expan¬ 
sion coefficients, and kb are definite wavenumbers spec¬ 
ified below. The prefactor is chosen so that for a single 
k mode {wi = l,Wb>i = 0) at small k (A:i?rms ‘C 1) the 
generalized eccentricity approaches £n,n 

lim e„{p(r)} = (13) 
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At small k, we expand the Jn{kr) and find 
— ^n,n 4“ Cl Sn^n+2^ 

where ci = —(fci?rnis/2)^/(l + n). Thus the functional 
form of p{r) adopted here yields a tunable linear com¬ 
bination of the eccentricities in Eq. (9), but the wave 
number parameter regulates the behavior at large r. Fur¬ 
ther motivation and discussion of this basis set for p(r) 
is given in our previous work [10]. 

We have found that an approximately optimal radial 
weight can be found by using only two well chosen kh 
values for the Fourier expansion in Eq. (12). Including 
additional k modes in the functional form of p{r) does 
not significantly improve the predictive power of the gen¬ 
eralized geometric eccentricity. For the two k modes we 
required (somewhat arbitrarily) that the ratio of k values 
would be fixed to the ratio of the first two Bessel zeros: 


With this choice our basis functions were orthogonal in 
the interval [0, i?o], where ki = jn,ilRo- We then ad¬ 
justed Ro to maximize the correlation coefficient between 
e„{p(r)} and the flow To account for changing sys¬ 
tem size with centrality, we used a fixed Ro/Rrms ratio. 
In most cases we used Ro/Rrms ~ 3.0, but for all directed 
flow components (^J^^ and ^P^) and the second elliptic 

flow component (Q found that Ro/Rrms ~ 2.0 op¬ 

timized the correlation between the flow and the geome¬ 
try. 

Ultimately, the assumption that the amplitude and 
phase of the flow is determined at least approximately by 
initial eccentricity, £„{p(r)}, is based on linear response. 
If nonlinear physics becomes important (as in the case 
of U 4 and us) then the predictors should be modified to 
incorporate this physics (see below and Ref. [20]). Thus, 
below we will refer to the £„{p(r)} (with an optimized 
radial weight) as the best linear predictor and incorpo¬ 
rate quadratic nonlinear corrections to the predictor as 
needed. 


III. RESULTS 
A. Radial flow 

Radial flow (or Vo(pt)) is the first term in the Fourier 
series and is by far the largest harmonic. Traditionally, 
the experimental and theoretical study of the fluctuations 
of Vb(PT) (i-e., multiplicity andpr fluctuations) has been 
distinct from elliptic and triangular flow. There is no 
reason for this distinction. 

Examining the scaled Vq(jpt) eigenvalues shown in 
Fig. 2(a), we see that there are two large principal com¬ 
ponents. The first principal component is sourced by 
multiplicity fluctuations, i.e., the magnitude of Vq(pt) 


fluctuates (but not its shape) due to the impact parame¬ 
ter variance in a given centrality bin. Corroborating this 
inference. Fig. 2(a) shows the momentum dependence of 
the leading principal component, which is approximately 
flat.^ Clearly this principal component is not particu¬ 
larly interesting, and the PCA procedure gives a practical 
method for isolating these trivial geometric fluctuations 
in the data set. The second principal component is of 
much greater interest, and shows a linear rise with pT 
that is indicative of the fluctuations in the radial flow 
velocity of the fluid [9] . 

In early insightful papers [23, 24], the fluctuations in 
the flow velocity (or mean px) were associated with the 
fluctuations in the initial fireball radius. These radial 
fluctuations are well described by both the eccentrici¬ 
ties (eq.o, £ 0 , 2 ), Eq. (9), and the optimized eccentricity 
eo{p{r)}, Eq. (II). Therefore, as seen in Fig. 2(b), the 
subleading flow signal is strongly correlated with these 
linear geometric predictors. 

Also shown in Fig. 2(b) is the correlation between sub- 
f 2) 

leading radial flow ' and mean transverse momentum 
fluctuations around the average 

SpT = [pt] - (br]) • (16) 

Indeed, the subleading radial flow correlates very well 
with mean momentum fluctuations in all centrality bins. 

In the next sections we will study the nonlinear mixing 
( 2 ) 

between the radial flow Q ' and all other harmonics. 

B. Elliptic flow 

1. Nonlinear mixing and elliptic flow 

We next study the fluctuations of V 2 {pt) as function 
of centrality. As seen in Fig. 3, the principal component 
spectrum of elliptic flow in central collisions consists of 
two nearly degenerate subleading components in addi¬ 
tion to the dominant leading component. This degen¬ 
eracy is lifted in more peripheral bins. Comparing the 
Pt dependence of the principal flows shown in Figs. 4(a) 
and 4(b), we see that going from central (0-5%) to pe¬ 
ripheral (45-50%) collisions, the magnitude of the sec¬ 
ond principal component increases in size and its momen¬ 
tum dependence changes dramatically. By contrast, the 
growth of the third principal component is much more 
mild. This strongly suggests that the average elliptic ge¬ 
ometry is more important for the subleading than the 
subsub-leading mode. 


^ There is a small upward tending slope in our simulations of this 
component, because multiplicity and mean px fluctuations only 
approximately factorize into leading and subleading principal 
components. Using different definitions of centrality bins could 
perhaps make this separation cleaner. 
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FIG. 2. (a) The pr dependence of the principal components of radial flow normalized by the average multiplicity, = 

Vq“^(Pt)/ {dN/dpr) ■ (b) The Pearson correlation coefficient [Eq. (8)] between the subleading radial flow and various predictors 
versus centrality. The best linear predictor is described in Sect. IIC. 
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FIG. 3. The magnitudes of the principal components of ellip¬ 
tic flow, 11 ^ 2 “^ II, versus centrality [see Eq. (6)]. 


To find a geometrical predictor for the sub- and sub- 
sub-leading modes we first tried the best linear predictor 
In Fig. 5(a) (the red circles), we see that the 
correlation coefficient between this optimal linear predic¬ 
tor and the subleading flow signal drops precipitously as 
a function of centrality. As we will explain now, this 
is because nonlinear mixing becomes important for the 
subleading mode. 

The ellipticity of the almond shaped geometry in pe¬ 
ripheral collisions is traditionally parametrized by eccen¬ 
tricity £ 2,2 and it serves as an excellent predictor for 


the leading elliptic flow. However, £ 2,2 does not com¬ 
pletely fix the initial geometry, and the radial size of the 
fireball can fluctuate at fixed eccentricity. As explained 
in Sect. Ill A, the radial size fluctuations modulate the 
momentum spectrum of the produced particles, and for 
a background geometry with large constant eccentricity 
this generates fluctuations in the pp dependence of the 
elliptic flow, i.e., subleading elliptic flow. This sublead¬ 
ing flow lies in the reaction plane following the average 
elliptic flow, but its sign (which is determined by Spr) is 
uncorrelated with £ 2 , 2 - 

The orientation of the reaction plane in peripheral bins 
is strongly correlated with the integrated V 2 or the lead¬ 
ing elliptic principal component , while the mean pp 
fluctuations are tracked by the subleading radial flow 
component Therefore we correlated the sub- and 

sub-sub-leading elliptic flows with the product of the 
leading elliptic and radial flows, i.e. we computed the 
correlation coefficient in Eq. (8) with C 2 %ed = 
Examining Fig. 5(a) (the black line), we see see that the 
correlation between the subleading elliptic flow and the 
nonlinear mixing rises with centrality, as the correlation 
with best linear predictor drops. Examining Fig. 5(b) 
on the other hand, we see that the subsub-leading ellip¬ 
tic flow has stronger correlation with the initial geom¬ 
etry than the nonlinear mixing. Combining best linear 
geometric predictor and quadratic mixing terms in the 
predictor, i.e. 


dpled = <^2{p{r)} + ciC2^^Co^\ (17) 

we achieve consistently high correlations for all centrali¬ 
ties [the blue diamonds in Fig. 5(a) and (b)]. 
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FIG. 4. The pr dependence of the principal components of elliptic flow normalized by the average multiplicity, v^\pt) = 
V 2 °'\pt)/ {dN/dpr), for central (0-5%) and peripheral collisions (45-50%). 
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FIG. 5. Pearson correlation coefficient between the subleading elliptic flows and the best linear predictor [Eq. (11)] with and 
without the nonlinear mixing between the radial and leading elliptic flows, (^) (b) show the correlation coefficient 

for V 2 subleading and V 2 subsub-leading flows respectively. 


2. Dependence on viscosity 


Before leaving this section we will briefly comment 
on the viscosity dependence of these results. Figure 6 
shows a typical result for a slightly larger shear viscosity, 
p/s = 0.16. As discussed above, the subleading ellip¬ 
tic flow [i.e., the event-by-event fluctuations in V 2 (pt)] 
is a result of the linear response to the first radial exci¬ 
tation of the elliptic eccentricity, and a nonlinear mix¬ 


ing of radial flow fluctuations and the leading elliptic 
flow. In Fig. 6 we see that a slightly larger shear vis¬ 
cosity tends to preferentially damp the linear response 
leaving a stronger nonlinear signal. This is because the 
initial geometry driving the linear response has a signif¬ 
icantly larger gradients due to the combined azimuthal 
and radial variations. Thus in Fig. 6 the linear response 
dominates the subleading flow only in very central colli¬ 
sions. These trends with centrality are qualitatively fa¬ 
miliar from previous analyses of the effect of shear vis- 
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FIG. 6. Pearson correlation coefficients for the subleading 
elliptic flow at viscosity over entropy ratio r]!s = 0.16. Dashed 
lines repeat r\ls = 0.08 results from Fig. 5(a) for the ease of 
comparison. 

cosity on the nonlinear mixing of harmonics [21, 25]. 

C. Triangular and directed flows 

Triangular flow was extensively studied in our previous 
work [10]. For the sake of completeness we relegate sev¬ 
eral comparative plots to the Appendix. Previously, we 
constructed an optimal linear predictor e 3 {/o(r’)} for the 
subleading triangular mode which characterizes the radi¬ 
ally excited triangular geometry. As shown in Fig. 7(b), 
adding the nonlinear mixing term to the best lin¬ 

ear predictor marginally improves the already good cor¬ 
relation with the subleading flow in peripheral collisions. 

Directed flow exhibits many similarities to triangular 
flow. Specifically the subleading directed flow is rea¬ 
sonably well correlated with the optimal linear predic¬ 
tor, characterizing the radially excited dipolar geometry. 
Nonlinear mixing between the leading directed flow and 
the radial flow is unimportant [see Fig. 7(a)]. 

D. The n = 4 and n = 5 harmonic flows 

It is well known that the leading components of the 
n = 4 and n = 5 harmonics are determined by the non¬ 
linear mixing of lower order harmonics in peripheral col¬ 
lisions [5, 18-21]. 

For comparison with other works [20, 26], in the Ap¬ 
pendix in Figs. 13(d) and 14(d) we construct a predictor 
based on a linear combination of the eccentricities 

£ 4,4 + Cie 2 , 2 e 2,2 for 71 = 4, (18) 

£5,5 + Ci£2,2e3,3 for n = 5, (19) 


where here and below the coefficient ci is adjusted to 
maximize the correlation with the flow. This predictor is 
compared to a linear combination of the optimal eccen¬ 
tricity £„{p(r)} and the corresponding nonlinear mixings 
of the leading principal components 

£4 {p(r)} + for 71 = 4, (20a) 

£5 {/o(£)} + for n = 5. (20b) 

Both sets of predictors perform reasonably well, though 
the second set has a somewhat stronger correlation with 
the flow. 

Returning to the subleading components, we first cor¬ 
related with the best linear predictors, £ 4 {p(r)} and 
£ 5 {p(r)}, with the corresponding subleading flow signals. 
As seen in Fig. 8 (the red circles), the correlation de¬ 
creases rapidly with centrality, especially for V4. Moti¬ 
vated by Eq. (20) which predicts the event-by-event lead¬ 
ing V 4 and Us in terms V 2 and U 3 , we construct a predictor 
for the subleading V4 and U 5 in terms of the fluctuations 
of V2 and U 3 (see Secs. IIIBl and IIIC, respectively). 
The full predictor reads 

£4 {/o(r)} -b > for 77 = 4, (21a) 

£5 {p(?')} + + C 24 ^^C 2 ^^ for 77 = 5. (21b) 

Including the mixings between the subleading V2 and U 3 
and the corresponding leading components greatly im¬ 
proves the correlation in mid-central bins (the blue dia¬ 
monds). Finally, in an effort to improve the V4 predictor 
in the most peripheral bins we have added additional 
nonlinear mixings between the radial flow and the lead¬ 
ing principal components 

£4 {p{‘r)} + 01^2^^'' + for r7 = 4, 

(22a) 

£5 {p{r)} + Ci? 2 ^^Cr^ + C 2 ? 3 ^^d^^ + for n = h. 

(22b) 

As seen in Fig. 8 (a) (the grey line) the coupling to the 
radial flow improves the correlation between the sublead¬ 
ing V4 and the predictor in peripheral collisions. On the 
other hand, for U 5 , Fig. 8 (b), all of the information about 
the coupling to the radial flow is already included in 
Eq. (21b) and adding vq does not improve the correla¬ 
tion. 


IV. DISCUSSION 

In this paper we classified the event-by-event fluctu¬ 
ations of the momentum dependent Eourier harmonics 
Vnipr) for 77 = 0-5 by performing a principal component 
analysis of the two-particle correlation matrix in hydro- 
dynamic simulations of heavy ion collisions. The leading 
principal component for each harmonic is very strongly 
correlated with the integrated flow, and therefore this 
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FIG. 7. Pearson correlation coefficient between the subleading (a) directed and (b) triangular flows and the best linear 
predictor with and without radial flow mixing. 
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component is essentially the familiar Vnipr) measured in 
the event plane. The subleading components describe ad¬ 
ditional pt dependent fluctuations of the magnitude and 
phase of u„(pt)- This paper focuses on the physical ori¬ 
gins of the subleading flows, which are the largest source 
of factorization breaking in hydrodynamics. 

Our systematic study started by placing radial flow 
(the n = 0 harmonic) in the same framework as the 
other harmonic flows. We identified the subleading n = 0 
principal component with mean px fluctuations and con¬ 
firmed (as is well known [23, 24]) that these fluctuations 


are predicted by the variance of the radial size of the 
fireball. 

The subleading directed and triangular flows were 
shown to be a linear response to the radial excitations 
of the corresponding eccentricity of the initial geometry. 
In these cases a generalized eccentricity e„{p(r)} with 
an optimized radial weight (describing the radial excita¬ 
tion) provides a good predictor for the subleading flows 
(Fig. 7). This extends our previous analysis of V 3 to 
vi [ 10 ]. 

Next, we investigated the nature of the subleading el- 
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liptic flows. The principal component analysis reveals 
that in central collisions there are two comparable sources 
of subleading elliptic flow, but they have strikingly dif¬ 
ferent centrality dependence (see Figs. 3 and 4). In 
mid-peripheral collisions the first subleading component 
mainly reflects a nonlinear mixing between elliptic and 
radial flows, and this component is only weakly corre¬ 
lated with the radially excitations of the elliptic geome¬ 
try. The second subleading component in this centrality 
range is substantially smaller and more closely reflects 
the radial excitations. In more central collisions, how¬ 
ever, the nonlinear mixing with the average elliptic flow 
becomes small, and the sub and subsub-leading principal 
components become comparable in size. Thus, the rapid 
centrality dependence of factorization breaking in V 2 is 
the result of an interplay between the linear response to 
the fluctuating elliptic geometry, and the nonlinear mix¬ 
ing of the radial and average elliptic flows. 

This nonlinear mixing can be confirmed experimen¬ 
tally by measuring the correlations between the princi¬ 
pal components ^^ 2 ^^which is predicted in 
Fig. 5. The prediction is that three point correlation be¬ 
tween the subleading elliptic event plane, the mean px 
fluctuations, and the leading elliptic event plane dehned 
by the Q 2 vector, i.e., 

U^2^5ptQ^S 

/ J=, (23) 


changes rapidly from central to midperipheral collisions. 
This correlation is analogous to the three plane correla¬ 
tions such as {V^{V 2 V^)*) measured previously [5]. 

Finally, we studied factorization breaking in W 4 and 
U 5 . With the comprehensive understanding of the fluc¬ 
tuations of V 2 and U 3 described above, the corresponding 
fluctuations in U 4 and U 5 were naturally explained as the 
nonlinear mixing of subleading V 2 and V 3 with their lead¬ 
ing counterparts, together with linear response to the 
quadrangular and pentagonal geometries (see Fig. 8 ). 

The study of the fluctuations in the harmonics spec¬ 
trum presented here shows the power of the principal 
component method in elucidating the physics which drive 
the event-by-event flow. We hope that this motivates 
a comprehensive experimental program measuring the 
principal components and their correlations for n = 0 — 5. 
Such an analysis would clarify the initial state in typical 
and ultra-central events with unprecedented precision, 
and would strongly constrain the dynamical response of 
the quark gluon plasma. 

ACKNOWLEDGMENTS 

We thank Jean-Yves Ollitrault, Wei Li, and Jiangyong 
Jia for continued interest. We especially thank Soumya 
Mohapatra for simulating hydro events. This work was 
supported by the U.S. Department of Energy under Con¬ 
tract No. DE-FG-88ER40388 


[1] Ulrich Heinz and Raimond Snellings, “Collective flow and 
viscosity in relativistic heavy-ion collisions,” Ann. Rev. 
Nucl. Part. Sci. 63, 123-151 (2013). 

[2] Matthew Luzum and Hannah Petersen, “Initial State 
Fluctuations and Final State Correlations in Relativis¬ 
tic Heavy-Ion Collisions,” J. Phys. G41, 063102 (2014), 
arXiv;1312.5503 [nucl-th]. 

[3] Vardan Khachatryan et al. (CMS), “Evidence for trans¬ 
verse momentum and pseudorapidity dependent event 
plane fluctuations in PbPb and pPb collisions,” Phys. 
Rev. C92, 034911 (2015), arXiv: 1503.01692 [nucl-ex]. 

[4] Georges Aad et al. (ATLAS), “Measurement of event- 
plane correlations in ^sjvjv = 2.76 TeV lead-lead colli¬ 
sions with the ATLAS detector,” Phys. Rev. C90, 024905 
(2014), arXiv:1403.0489 [hep-ex]. 

[5] Georges Aad et al. (ATLAS), “Measurement of the 
correlation between flow harmonics of different order 
in lead-lead collisions at y^SiCv=2.76 TeV with the 
ATLAS detector,” Phys. Rev. C92, 034903 (2015), 
arXiv:1504.01289 [hep-ex]. 

[6] Jaroslav Adam et al. (ALICE), “Event shape engineering 
for inclusive spectra and elliptic flow in Pb-Pb collisions 
at y'sNN = 2.76 TeV,” arXiv: 1507.06194 [nucl-ex]. 

[7] L. Adamczyk et al. (STAR), “Long-range pseudo¬ 
rapidity dihadron correlations in d-|-Au collisions at 

= 200 GeV,” Phys. Lett. B747, 265-271 (2015), 
arXiv;1502.07652 [nucl-ex]. 


[8] A. Adare et al. (PHENIX), “Measurements of elliptic and 
triangular flow in high-multiplicity ®He-|-Au collisions at 

= 200 GeV,” Phys. Rev. Lett. 115, 142301 (2015), 
arXiv: 1507.06273 [nucl-ex]. 

[9] Rajeev S. Bhalerao, Jean-Yves Ollitrault, Subrata Pal, 
and Derek Teaney, “Principal component analysis of 
event-by-event fluctuations,” Phys. Rev. Lett. 114, 
152301 (2015), arXiv:1410.7739 [nucl-th]. 

[10] Aleksas Mazeliauskas and Derek Teaney, “Sublead¬ 
ing harmonic flows in hydrodynamic simulations of 
heavy ion collisions,” Phys. Rev. C91, 044902 (2015), 
arXiv:1501.03138 [nucl-th]. 

[11] Fernando G. Gardim, Frederique Grassi, Matthew 
Luzum, and Jean-Yves Ollitrault, “Breaking of factor¬ 
ization of two-particle correlations in hydrodynamics,” 
Phys. Rev. C87, 031901 (2013), arXiv:1211.0989 [nucl- 
th]. 

[12] K. Aamodt et al. (ALIGE), “Harmonic decomposition of 
two-particle angular correlations in Pb-Pb collisions at 
yslEV = 2.76 TeV,” Phys. Lett. B708, 249-264 (2012), 
arXiv: 1109.2501 [nucl-ex]. 

[13] Georges Aad et al. (ATLAS), “Measurement of the 
azimuthal anisotropy for charged particle production 
in ^/snn = 2.76 TeV lead-lead collisions with the 
ATLAS detector,” Phys. Rev. C86, 014907 (2012), 
arXiv: 1203.3087 [hep-ex]. 

[14] Georges Aad et al. (ATLAS), “Measurement of long- 









10 


range pseudorapidity correlations and azimuthal har¬ 
monics in ^snn = 5.02 TeV proton-lead collisions with 
the ATLAS detector,” Phys. Rev. C90, 044906 (2014), 
arXiv;1409.1792 [hep-ex]. 

[15] Ulrich Heinz, Zhi Qiu, and Chun Shen, “Fluctuating flow 
angles and anisotropic flow measurements,” Phys. Rev. 
C87, 034913 (2013), arXiv: 1302.3535 [nucl-th]. 

[16] Igor Kozlov, Matthew Luzum, Gabriel Denicol, Sangyong 
Jeon, and Charles Gale, “Transverse momentum struc¬ 
ture of pair correlations as a signature of collective behav¬ 
ior in small collision systems,” arXiv:1405.3976 [nucl-th]. 

[17] B. Alver, M. Baker, C. Loizides, and P. Steinberg, “The 
PHOBOS Glauber Monte Carlo,” arXiv:0805.4411 [nucl- 
ex]. 

[18] Nicolas Borghini and Jean-Yves Ollitrault, “Momentum 
spectra, anisotropic flow, and ideal fluids,” Phys. Lett. 
B642, 227-231 (2006), arXiv:nucl-th/0506045 [nucl-th]. 

[19] Zhi Qiu and Ulrich W. Heinz, “Event-by-event shape and 
flow fluctuations of relativistic heavy-ion collision fire¬ 
balls,” Phys. Rev. C84, 024911 (2011), arXiv: 1104.0650 
[nucl-th] . 

[20] Fernando G. Gardim, Frederique Grassi, Matthew 
Luzum, and Jean-Yves Ollitrault, “Mapping the 
hydrodynamic response to the initial geometry in 
heavy-ion collisions,” Phys. Rev. C85, 024908 (2012), 
arXiv:1111.6538 [nucl-th]. 

[21] Derek Teaney and Li Yan, “Non linearities in the har¬ 
monic spectrum of heavy ion collisions with ideal and 
viscous hydrodynamics,” Phys. Rev. C86, 044908 (2012), 
arXiv:1206.1905 [nucl-th]. 

[22] Mikko Laine and York Schroder, “Quark mass thresh¬ 
olds in QCD thermodynamics,” Phys. Rev. D73, 085009 
(2006), arXiv;hep-ph/0603048 [hep-ph]. 

[23] Wojciech Broniowski, Mikolaj Chojnacki, and Lukasz 
Obara, “Size fluctuations of the initial source and the 
event-by-event transverse momentum fluctuations in rel¬ 
ativistic heavy-ion collisions,” Phys. Rev. C80, 051902 


(2009), arXiv:0907.3216 [nucl-th]. 

[24] Piotr Bozek and Wojciech Broniowski, “Transverse- 
momentum fluctuations in relativistic heavy-ion colli¬ 
sions from event-by-event viscous hydrodynamics,” Phys. 
Rev. C85, 044910 (2012), arXiv:1203.1810 [nucl-th]. 

[25] Zhi Qiu and Ulrich Heinz, “Hydrodynamic event-plane 
correlations in Pb-|-Pb collisions at ^/s = 2.76ATeV,” 
Phys. Lett. B717, 261-265 (2012), arXiv: 1208.1200 
[nucl-th] . 

[26] Fernando G. Gardim, Jacquelyn Noronha-Hostler, 
Matthew Luzum, and Frdrique Grassi, “Effects of viscos¬ 
ity on the mapping of initial to final state in heavy ion col¬ 
lisions,” Phys. Rev. C91, 034902 (2015), arXiv:1411.2574 
[nucl-th] . 


APPENDIX: LIST OF FIGURES 

Here we present a comprehensive catalog of plots for 
each harmonic n = 0-5. Centrality dependence of flow 
magnitudes for n = 0, appearing as Fig. 3 in the text 
above, is repeated here as Fig. 9(a), and analogous plots 
for other harmonics are given in Figs. 10-14(a). The pt 
dependence of normalized principal components for ra¬ 
dial and elliptic flows in central (0-5%) collisions shown in 
Figs. 2(a) and 4(a) are reproduced as Figs. 9(c) and 11(c) 
and complemented with Figs. 10(c) and 12-14(c). Addi¬ 
tionally, Figs. 9-14(b) depict the same principal compo¬ 
nents, but without normalization by average multiplicity 
{dN/dpT)- Finally, in the paper we showed the Pearson 
correlation coefficients for the subleading flows for each 
harmonic n = 0-5 in Figs. 2(b),5(a),7(a), 7(b), 8(a) and 
8 (b), while in this appendix we show results for both lead¬ 
ing and subleading flows in the series of figures Figs. 9- 
14(d) and Figs. 9-14(e). 
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